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Abstract Wc present a new improved version of our force-free electrodynamics 
(FFE) numerical code in spherical coordinates that extrapolates the magnetic 
field in the inner solar corona from a photospheric vector magnctogram. The 
code satisfies the photospheric boundary condition and the condition V • B = 
to machine accuracy. The performance of our method is evaluated with standard 
convergence parameters, and is found to be comparable to that of other nonlinear 
force-free extrapolations. 
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1. Introduction 

Force-free electrodynamics (hereafter FFE) is a formal name for time-dependent 
electromagnetism in an ideal plasma with negligible inertia and gas pressure 
(/? s=s 0). The formalism of FFE has been developed for various relativistic as- 
trophysical applications where the plasma supports electric currents and electric 
charges (pulsars, astrophysical jets, gamma-ray bursts, etc.). The equations of 
FFE are Maxwell's equations with nonzero electric currents, 

$E <9B 

— = cV x B - 4ttJ , — = -cV x E , (1) 
at at 

complemented by the divergence-free, ideal MHD, and force-free conditions 

V • B = , E • B = , /9 C E + - J x B = . (2) 

c 

Here, J and p c = (47r) _1 V • E are the electric current and charge densities, 
respectively. One can solve the above set of equations for J and thus express the 
electric current density as a function of the electric and magnetic fields, namely 

T c v F !l2iJi 4. c (B-VxB-E-VxE) 
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( jGruzinov, 1999[ ). One can then numerically integrate Maxwell's equations to 
follow the temporal evolution of any force-free ideal-MHD system from an ini- 
tial to a final configuration. In the particular case when E *— ? 0, the final 
configuration is (in general) a non- linear force-free magnetostatic field, i.e. it 
satisfies the relations 

J = -£-V xB, V-B = 0, JxB = 0. (4) 

47T 

The FFE method has been applied successfully in two main astrophysical ap- 



plications; pulsars ( Spitkovsky 2006; Kalapotharakos and Contopoulos I2009[l . 
and the solar corona fContopoulos, Kalapotharakos, and Georgoulis 120111 here- 



after CKG). The numerical codes developed for these problems are Cartesian 
finite-difference time-domain (FDTD) codes. They are based on the staggered 
mesh algorithm of lYeel (I1966[) where in every computational cell, electric field 
components arc defined in the middle of and along cell edges, whereas magnetic 
field components arc defined in the middle of and perpendicular to cell faces. It is 
very important to notice that this prescription satisfies the divergence-free con- 
dition V • B = to machine accuracy |Kalapotharakos and Cont opoulos (2009) 
introduced also an outer absorbing non-reflecting boundary condition in the 
form of a perfectly matched layer (PML). CKG implemented the FFE method 
in the study of the inner solar corona by evolving the photospheric boundary 
condition toward a given distribution for the radial component of the magnetic 
field. This is achieved with the introduction of photospheric electric fields which 
gradually die out as the required photospheric condition is approached. During 
that evolution, the FFE equations arc numerically solved in the solar corona, 
and in the limit when electric fields die out everywhere, a non-linear force-free 
magnetostatic configuration emerges. 

There are several problems with this approach: a) The method is based on 
knowledge only of the radial component of the photospheric magnetic field; 
therefore the final solution is not unique. Different initial conditions and different 
photospheric evolutionary paths yield different magnetostatic configurations, b) 
Numerical dissipation cannot be avoided, and therefore, the longer a particular 
implementation takes to converge to a certain acceptable accuracy, the smaller 
the final remaining electric currents in the corona, and the closer the final solu- 
tion approaches the potential (current-free) equilibrium, c) The numerical grid is 
Cartesian. This works well around the poles, but is unsuitable to appropriately 
model the solar photosphere in general. 

For all of the above reasons, we decided to improve our code, and to incor- 
porate the photospheric boundary condition supplied from a full vector magnc- 
togram . 



2. Implementation of the Vector Magnetogram Boundary Condition 

The new improved version of our code is written in the natural coordinates for the 
solar corona, namely heliocentric spherical ones (r, 0,<fi). In order to satisfy the 
photospheric boundary condition exactly, we altogether skip the photospheric 
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relaxation algorithm described in CKG, and begin instead with the following 
'sea urchin'-likc initial magnetic field configuration: 

B r (r,9^;t = 0)=B r (9^)(r Q /r) 2 , 



B e (r 0) M;* = O) = B*(M) , 
B 4 {r Q ,e,4>;t = 0) =B^{6,4>) , and 

B e (r > r 0! 0, 0; t = 0) = B^{r > r , 9, 0; t = 0) = . 

This configuration is obviously divergence-free, and satisfies the photospheric 
boundary condition B(r ,(9, <j>) = B{6,4>) exactly. We also set E = on the 
photosphere at all times. Notice that the algorithm that we implement (Yce, 
1966) satisfies V • B = to machine accuracy only if the initial configura- 
tion is divergence-free. In other words, the algorithm does not implement the 
divergence-free condition. It inherits and preserves it. Notice also that the initial 
configuration is filled with electric currents in the 9- and ^-directions, but not 
in the radial one. 

We then proceed as in CKG with a numerical integration of Equations ([1]) 
and (j3J). Initially, the magnetic field configuration becomes torn through nu- 
merical reconnection. The sea urchin-like configuration gradually disappears. 
Transient electromagnetic-type waves travel through the integration volume and 
when they reach the outer radial boundary they are absorbed by the PML (see 
below). Eventually, this transient activity dies out, electric fields are observed 
to converge to E q, and a non- linear force- free magnetostatic equilibrium 
is gradually reached. Notice that during this evolution, the divergence-free and 
vector photospheric boundary conditions are everywhere satisfied. 

Our method can be implemented in both global and local magnetic field 
extrapolations. In global extrapolations, the numerical integration volume lies 
between an inner radius r m = r (the photosphere), and an outer radius r out that 
corresponds to the distance beyond which the solar wind becomes dynamically 
important. The integration volume corresponds to the lower solar corona where 
we assume that force-free magnetostatic conditions apply to a certain extent. 
In the test simulations presented in the next section, we took r ou t = 2r . In 
general, r out may vary between 2r and 3r , depending on the strength of the 
solar wind. Beyond that radius, we implement a PML non-reflecting absorbing 
layer of thickness O.5r . In the polar direction, the simulation is limited between 
a minimum and a maximum polar angle # m in = 20° and 9 max = 160°, in order 
to avoid the polar singularity of the spherical coordinate system. In realistic 
solar coronal magnetic field configurations, the polar regions are covered by 
coronal holes, and therefore, their exclusion does not significantly affect the ex- 
trapolation. Finally, in the azimuthal direction we implement periodic boundary 
conditions. The test simulations presented in the next section are all global 
extrapolations except for the last one that is limited in the azimuthal extent. 
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3. Evaluation of Numerical Code 

In order to evaluate the performance of our code, we implement vector photo- 
spheric boundary conditions derived from standard solutions that are used as 
benchmarks, and evaluate the convergence of the extrapolation based on stan- 
dard convergence parameters introduced in Wheatland, Sturrock, and Roumeliotisj 



(pOOO)) . |Schrijver et al.\ ([2006]) . and |Amari, Boulmezaoud, and Aly| (|2006|) fsee 



the Appendix for their detailed expressions) . The benchmark solutions are known 
either analytically (dipole), or quasi-analytically (|Low and Loul [19901 hereafter 
LL). In most previous implementations, benchmark boundary conditions on the 
xy-plane were obtained by displacing the LL solution some distance / below the 
sy-plane, and by rotating it by an angle $ with respect to the y-axis (see the 
original LL paper for details). To the best of our knowledge, the only previous 
work where the LL solutions were applied in a coronal magnetic field configura- 
tion over a curved photospheric boundary is |Tadesse, W icgclmann , ~and Inhester| 
(200SJ). We thus decided, for the sake of comparison, to also implement their 
boundary condition. We investigated the following test cases: 

• Case DF: a dipole magnetic field off-axis by 0.5r© 

• Case LL11: LL solution with n = 1, m = 1, 1 = 0.25r© (displacement below 
the cry-plane), $ = 7r/10 (as in |Tadesse, W icgclm ann, and Inhester| |2"0"0"9")) 

• Case LL13: LL solution with n = 1, m = 3, I = O.lr (displacement below 
the j/z-plane), $ = 4ir/5 

• Case LL31: LL solution with n = 3, m = 1, I = 0.1r Q (displacement below 
the yz-planc), $ = 4ir/5 

Case DF is current-free (potential), but the LL solutions are not. The radial 
resolution of these particular simulations is O.lr©, and the angular resolution 
4° x 4° (heliocentric). We set m i n = 20° and 6* max = 160°. Our computational 
grid in (r, 8, </>) has a resolution of 10 x 34 x 90. 

In Table 1, we plot the values of our convergence parameters in a subregion 
of our computational volume (we chose r Q < r < 1.5r , 30° < 9 < 60°, 
30° < <fi < 60°) in the force- free extrapolations obtained for the above four 
benchmark solutions. The first four parameters characterize the convergence 
to the benchmark solution (unity corresponds to a perfect match). The last 
two parameters characterize how well the extrapolation satisfies the force-free 
and divergence-free conditions respectively. In all cases, the divergence-free and 
boundary conditions are satisfied to machine accuracy. The force-free condition 
is also very important for us. We evolve the numerical simulation till we obtain a 
minimum of the parameter CWsin at a level roughly below 0.15. This guarantees 
that the solutions that we present are indeed non-linear force-free magnetostatic 
configurations that satisfy the given vector magnetogram boundary conditions. 
Notice that the force- free parameter (Equation (fTO)) ) does not apply in potential 
extrapolations where J — > 0. 

The performance of our code is satisfactory, which suggests that the present 
numerical method is promising. In some benchmark cases, it outperforms existing 
Cartesian grid extrapolations (see e.g. Jiang and Feng, 2012]), and is certainly on 



a par with the spherical grid extrapolation of Tadesse, Wicgclm ann, and Inhesterj 
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Table 1. Convergence parameters for various benchmark extrapolations on a computa- 
tional grid with a resolution in (r, 8, <f>) of 10 X 34 X 90. 



Case 


Cvec 


Ccs 


K 


E 'm 


CWsin 


/ 




Comments 


DF 


0.989 


0.986 


0.845 


0.821 




10" 


9 


Current-free 


LL11 


0.993 


0.987 


0.895 


0.873 


0.087 


io- 


8 


as in Tadesse et al. (2009) 


LL13 


0.915 


0.924 


0.403 


0.389 


0.109 


io- 


8 


Different solution 


LL31 


0.981 


0.978 


0.772 


0.715 


0.110 


10" 


7 


Non-global in azimuth 



(2009) (compare our case LL11 with Case 1 in their low resolution grid). In some 
cases, though, the resulting extrapolation does not converge to the benchmark 
solution in the volume under investigation, as manifested for example in case 
LL13 where two convergence parameters differ significantly from unity. We offer 
two possible explanations for this result: a) An extrapolation where a signif- 
icant amount of magnetic flux crosses the outer boundaries of the numerical 
integration region is influenced by the presence of those boundaries more than 
another extrapolation where most magnetic field lines are contained within the 
computational volume, b) A non-linear extrapolation is not unique when current 
sheets are allowed to develop in the solution. This needs to be confirmed with 
more detailed evolutionary full MHD extrapolations. Notice that our force-free 
code handles current sheets as tangential magnetic field discontinuities with 
\B\ continuous everywhere (see |Contopoulos, Kazanas, and Fendt| ()1999p and 



Kalapotharakos and Contopoulos (|2009|) for details). Finally, we noticed that in 



a few cases the extrapolation did not converge, probably due to an incompat- 
ibility between the photospheric vector boundary condition and the conditions 
at the outer boundaries of the computational volume. As an example, the global 
extrapolation failed in case LL31; however, when we limited the azimuthal 
integration region to 20° < <j) < 160° the extrapolation was successful. The 
convergence properties of our method certainly need further investigation. 

The resulting solutions can be visually compared with the corresponding 
benchmark solutions in Figure 1 where field lines originating from the same 
photospheric positions arc plotted. Wc notice that the extrapolation produces 
field lines that are visibly more stretched-out radially than the benchmark so- 
lution. This is particularly obvious in the case of a current-free displaced dipole 
and in the field lines that cross the outer radial boundary at r = 2r . This result 
is an artefact of our method that qualitatively mimics the effect of the solar wind 
beyond the outer radial distance. 



4. Prospects for the Future 

The implementation of the photospheric boundary condition from full vector 
magnctograms in the new version of our code improves the nonlinear force-free 
extrapolation of CKG in several respects. Firstly, it generates a more physical 
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solution (closer to reality). Secondly, it resolves the non- uniqueness problem of 
the original version where different initial conditions and different photosphcric 
evolutions yielded different magnetostatic configurations. And thirdly, it docs 
not allow the extrapolation to 'degenerate' to the potential configuration. Our 
test runs suggest that the method is promising, as it outperforms existing ex- 
trapolations in some benchmark cases. We therefore expect that it has significant 
potential in the study of the lower solar corona and of individual active regions. 

We are currently running our code on a desktop PC. Our global coronal mag- 
netic field extrapolations were obtained with a heliocentric angular resolution of 
4° x 4°, on an (r, 9, (f) spherical computational grid of 10 x 34 x 90, but this can 
certainly be improved. Our goal is to parallelize the code to run in a supercom- 
puter. In a future implementation we will also improve the performance of the 
PML outer boundary by fine-tuning its parameters (our current implementation 
allows for some minor amount of reflection of transient clcctromagnetic-type 
waves). 

An ambitious goal is to apply our method to investigate particular active 
regions with minimum influence from the simulation outer boundaries on the 
extrapolation. As we saw, the effect of the outer boundaries differs from case 
to case, and we can only give an empirical rule that the assumed volume must 
be about double the size of the particular volume under investigation in all 
directions (radial, polar, and azimuthal). Our goal for the future is to solve 
for the magnetic field of a particular region as part of (i.e. together with) a 
global coronal magnetic field extrapolation. This can be achieved in practice by 
introducing adaptive mesh refinement (AMR) around the region under investi- 
gation and around neighboring regions that may interact with that region. In 
practice, we will need about five adaption levels in order to match the available 
heliocentric angular resolution in particular active regions (about 0.05° x 0.05°) 
to the resolution that can be practically achieved on a single PC in a global 
magnetic field extrapolation (about 2° x 2°). Thus, the extrapolation for one 
particular active region will take into account the presence of neighboring active 
regions, as well as the global coronal magnetic field structure. 

One final goal is to introduce finite (non-infinite) conductivity a in the code 
by modifying the FFE expression for the electric current density as 

E x B 

J = /°c— ga~ + , (5) 

We plan to investigate various nonlinear prescriptions for a (e.g. weighted by the 
value of J) as we are currently doing in our study of the pulsar magnetosphere 



(see Kalapotharakos et al.\ ((2012') for details). Finite conductivity will allow for 



the presence of electric fields along the magnetic field, i.e. it will introduce 
volume energy dissipation J • E. We expect that this approach will give us direct 
hints about where and how magnetic field energy is released in the solar corona. 
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Appendix 

Convergence Parameters 



We evaluate the following parameters introduced in Wheatl and, Sturrock, and Roumeliotisj 



pTOD)) , |Schrijver etaJ] (|2005|l . and |Amari, Boulmezaoud, and Aly| ([20"0"o]) : the 
vector correlation C vcc 



i \ i i 



\ 1/2 
|2 I 



(6) 



the parameter Ccs 



and the normalized and mean vector error £7^ and E' m 

K^l-^IBi-bil/^lBil , (8) 

i i 

where and are the extrapolation and reference fields, respectively, i denotes 
grid points, and N is the total number of grid points in the sub-region of our 
simulation where we evaluate the convergence of the code. The force-free and 
divergence-free conditions are estimated using the following parameters: 

CWsin ^'V B f» B *' , (10) 

/ = (ii) 
J ~ N ^ 6|Bi|/d ' 1 ' 

where d is the minimum grid spacing. Those may be evaluated over the total 
integration volume, or over a smaller internal volume. 
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Figure 1. Comparison between benchmark solutions (left) and their corresponding extrapo- 
lations (right) for the four cases of Table 1. In case LL11, the disk center corresponds to 180° 
longitude. In the other three cases, it corresponds to 90° longitude. Dashed lines mark the 
simulation photospheric boundaries. 
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